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Abstract. We present a reduction pipeline for CCD {charge-coupled device) images which was built to search for 
variable sources in highly crowded fields like the M31 bulge and to handle extensive databases due to large time 
series. We describe all steps of the standard reduction in detail with emphasis on the realisation of per pixel error 
propagation: Bias correction, treatment of bad pixels, flatfielding, and filtering of cosmic rays. The problems of 
conservation of PSF {point spread function) and error propagation in our image alignment procedure as well as 
the detection algorithm for variable sources are disc ussed: We build differe nce images via image convolution with 
a technique called OIS {optimal image subtraction, Alard & Lupton, 1998), proceed with an automatic detection 
of variable sources in noise domina ted images and finally apply a PSF-fitting, relative photometry to the sources 
found. For the WeCAPP project (Riffeser et al., 2001) we achieve 3a detections for variable sources with an 
apparent brightness of e.g. m — 24.9 mag at their minimum and a variation of Am — 2.4 mag (or m = 21.9 mag 
brightness minimum and a variation of Am = 0.6 mag) on a background signal of 18.1 mag/arcsec'^ based on a 
500 s exposure with 1.5 arcsec seeing at a 1.2 m telescope. The complete per pixel error propagation allows us to 
give accurate errors for each measurement. 



Key words. Methods: data analysis - Methods: observational 
propagation - Techniques: optimal image subtraction 



Techniques: image processing - Techniques: error 



1. Introduction 

Astronomical imaging in optical wavebands is performed 
nearly exclusively with charge- coupled devices^ today. 
Despite the numerous advantages of modern CCDs, their 
images still have to be corrected for a couple of disturb- 
ing influences and effects before one can base advance in 
science on them. Here, we will focus on the problems aris- 
ing with optical, ground based imaging and time-series 



observiiiLlous Lu find and lueasLire variable sources either 
hidden 111 a bright background (e.g. a variable star hi its 
host galaxy) or a crowded field or even in a combination 
of both. 

The search for variable objects with common photom- 
etry methods becomes very ineffective in crowded fields 
because of blending. Phillips fc Davis (1995| ) show algo- 
rithms for registering, matching the point spread functions 
(PSFs), and matching the intensity scales of two or more 
images in order to detect transient events. Tomaney & 
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The history of CCDs in astronomy a n d a basic de s cription 
of them can be found in 
(1990D ; |Mackay {198^ . 



McLean (1997); Bull (1991); Jacoby 



Crotts (1996) propose a method called Difference Image 
Analysis (DIA) where the point spread function (PSF), 
describing the projection of a point source onto the im- 
age plane, is matched by calculating a convolution kernel 
in Fourier space. The application to real data has been 



shown for Galactic Microlensing (Alcock et al., 1999) as 
well as for microlensing in M31 ( Crotts et al., 1999| ). 
A new method for Optimal Image Subtraction (OIS) 
of two images has been designed by Alard fc Luptor] 
(1998). They derive an optimal kernel solution from a 



simple least-squares analysis using all pixels of both im- 
ages. This method has been used successfully by different 
projects ( OGLE,|Wozniak, 200C|; MOA, [Bond et al., 2001 
DIRECT, iMochejska et al., 20011; etc.). 



We have expanded OIS with a standard reduction 
pipeline, a finding algorithm for variable sources, a PSF- 
fitting photometry, and per pixel error propagation for all 
steps of image processing. We are able to give accurate 
errors for each photometric measurement of a variable 
source beyond an estimate based on the noise in a sin- 
gle difference image. Our software has been developed to 
deal with time series obtained with the 0.8 m Wendelstein 
telescope and the Calar Alto 1.23 m telescope to monitor 
variable stars and find microlensing events towards the 
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origin property 

Detector: Photosensitive area and additional "borders" (prescan, postscan, overscan), geometric variation of pixel 

size and edge pixels, pixel sensitivity variations and pixel defects (cold pixel, hot pixel, trap), sub-pixel 
quantum efficiency variations, charge transfer efficiency (CTE), linearity range and saturation level. 

Electronics: Bias, gain, ADC width, sampling (of charges), thermal noise. 

Instrument: Dust in optics and optic distortions, optical scale and detector pixel size (and their ratio to the typical 

seeing - spatial sampling), file format (e.g. FITS ^), and information beyond raw image (header keywords). 

Environment: Signals originating from particle events (cosmic rays), varying meteorological observing conditions (seeing), 
and other atmospheric effects like sky illumination (moon) and extinction. 

Table 1. Properties of raw CCD images and their origin. 



M31 bulge in the first instance. However, we stress that 
many of the procedures presented here may also be applied 
directly to other types of astronomical imaging. 

The first part of this paper motivates the effort to im- 
plement per pixel error propagation. The second part gives 
a detailed description of our standard reduction for CCD 
frames. In the third part we present our image alignment 
procedure. The fourth part describes the image convolu- 
tion with OIS, the detection procedure for variable sources 
and the relative photometry of those sources. In the last 
section we give the results of performance tests on simu- 
lated images. 



2. Motivation for per pixel error propagation 

2.1. Properties of raw CCD images 

We have to consider all properties of raw CCD images 
(Tab. ^ before we can establish a reliable difference im- 
age analysis. When estimating photometric errors all these 
effects and their errors have to be considered again in ad- 
dition to the photon noise induced by the objects imaged. 
To account for those additional photometric errors the er- 
ror has to be calculated and propagated for each pixel and 
every data reduction step applied. 



2.2. Error propagation - an analytic example 

We show that the errors resulting e.g. from flatfield cal- 
ibration (Sect. 3.7) may be dominant for bright objects 
or inadequate flatfields which cannot always be avoided. 
Statistics for this effect on simulated images are shown in 
Sect. 3.1. (See Tab. |^ for common components of formulae 
used throughout this article.) 

One can estimate this effect assuming nj images, each 
with I{x,y) detected photons per pixel, and 5i(x,y) = 
\J I{x, y) error per pixel, so the preliminary error estimate 
for a per pixel added stack 



ni 



l{x,y) 



■[F{x,y) F{x,y) 



(1) 



Greisen et al., 198f 



Flexible Image Transport System, see Wells et al., I98I 



Grosbol et al., 1988 



Ponz et al., 1994 ^^nd NOST 100-2.0. 



Harten et al., 1988 



{x, y) = pixel coordinates, 
I{x,y) = value of pixel {x,y) in image /, 
Si{x,y) = corresponding absolute error, 
1,1,1 = sequence of indicators, that some reduction step 
has been applied to /, 
nj = integer number of e.g. images of type /, 
a = the root mean square of a sample. 

Table 2. Common components for the notation of formu- 
lae. 



(where is due to the fact that the Ij will vary, and 
F{x,y) « 1) is 



Sj{x,y) « Si{x,y) y/nj . 



(2) 



Including as further assumptions a flatfield error 
5f{x, y), a ratio of the relative errors between the flatfield 
and a single image given by 



S.{x,y) 



Siix,y) j &F{x,y) 
Iix,y) I F{x,y) 



(3) 



and Gaussian error propagation, we can calculate the im- 
pact of flatfleld calibration on the error Jj(a;, y) of the 
stack. The propagated per pixel error depends on the indi- 
vidual signal-to-noise ratios of the images and the flatfield, 
and on dithering. When stacking spatially undithered im- 
ages, the flatfield error SF{x,y) of the same pixel {x,y) 
adds to each individual image I{x, y), so the flatfleld error 
part of the stack is not an independent error. It is like 
stacking first and flatflelding subsequently: 



5j{x, y) w 5i{x, y) xjni [1 + 



ni 



(4) 



With spatially dithered (and digitally realigned) images 
we get (by neglecting the effects of the alignment proce- 
dure) 



5j{x,y) w 5i{x,y) W ti/ 1 + 



1 

e 



(5) 



because different and therefore independent flatfleld pixels 
add to the per pixel error of the stack. 

An example: We consider an extended and bright ob- 
ject (e.g. the M31 bulge), and twilight flatflelds. Here we 
get 1 < ^ < 3 because of the difficulties in getting twilight 
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flats discussed in Sect. 3.7. In a dithered, nj — 5 stack 



this yields an up to 41% error increment and even a 145% 
increment in an undithered, nj = 5 stack, both compared 
to the simple estimate of Eq. |^, i.e. with respect to the 
Poissonian noise in the images alone. 

Now we treat a complementary case: When imaging 
faint objects in empty fields, the lack of proper Hatfield^ 
pushes observers to build flatfields using the skylight in the 
(dithered) images. With I{x,y) = object(a;,y) + sky{x,y) 
and Hp of 71/ images usable for building the flat (np < nj 
because of the objects and cosmics in the field) we get 



object(x, y) 
sky(a:, y) 



for the error in a stack, using 



1 



5F{x,y) ^ 

F{x,y) y/npsky 



in Eq. 5. 



(6) 



(7) 



For Hp — 5 images, which can be used in every pixel for 
flatfielding, and with an object which is much fainter than 
the sky, this yields a 10% error increment (here, and below 
again, compared to the naive, Poissonian noise estimate 
of Eq. 1^). Assuming an object with flux comparable to the 
sky, the error increment already exceeds 18%. 

The per pixel propagation of errors gets even more 
important when performing multi-pixel approximations, 
as we do in some parts of our data reduction pipeline 
for the difference image analysis. When performing e.g. 
PSF-fitting one can enhance the accuracy by including 
the appropriate error weights as stated in Sect. 6.6. 



3. Standard reduction pipeline for images 

We start with the standard red uction for individual object 
frames. All effects of Sect. 2.1, which can be corrected in 



single images, will be discussed. We show a way how to 
compensate for these effects, and, in addition, how these 
compensations affect the total error budget of each pixel. 

3.1. Saturation 

Saturated pixels (due to ADC saturation or due to pixel 
charges at full well capacity) have to be marked as wrong 
value pixels. With a reasonable bias adjustment no pixel 
should have a zero value, so we set saturated pixels to 
zero as a first step. Since one cannot know if a saturated 
pixel has or has not bloomed, the four directly adjacent 
neighbouring pixels should also be treated like saturated. 
Depending on the specific properties of a CCD, the bloom- 
ing assumption may be restricted to the two pixels in-and- 
against the row-shift direction. As we perform the bias 
correction (Sect. |3.2| ) and the initial calculation of an er- 
ror frame (Sect. |3.3[ ) in one step, the zero marked pixels 
will still be recognised after the bias has been subtracted. 



It is nearly impossible to get twili ght flatfields for more 
than three filters in one night; see Sect. 3.7. 



Pixels with values beyond the linearity range of the CCD 
have to be corrected, or, if not possible, also have to be 
treated like saturated pixels. Since we do not have to treat 
images with values in the non linear regime we will not 
provide error propagation for this case. 

3.2. Bias correction 

The bias originates from the offset voltage of the CCD 
ADC. Depending on its properties we remove it following 
two subsequent approaches: 

— If the patterns in the bias frame are varying on short 
time scales, or there is no pattern at all (just ther- 
mal noise) , we subtract only the Kcr-clipped mean of a 
suitable part of the overscan (i.e. a part of the over- 
scan, which is identical in exposed frames and dark, 
not exposed frames) . Since the bias is an additive con- 
stant and mostly a small number the Kcr-clipped mean 
(k w 6, to get rid of cosmics) will be more accurate 
than the median. 

— If the bias patterns are reproducible, we subtract in 
addition a Kcr-median-clipped mean image of overscan 
corrected bias frames. (We check the reproducibility by 
looking for patterns in such a mean image; see Sect. 3.7 
for details on kct- median-clipping.) 

Three types of errors have to be considered when correct- 
ing for the bias: 

— Statistical error: The readout noise consists of the ther- 
mal noise of the amplifier, and, if present, of the dark 
current. 

— Systematic error: Unreproducible bias patterns may 
result from bad CCD electronics or insufficient elec- 
tronic shielding. (Strong radio immission may even 
penetrate an excellent shielding.) 

— Numerical error: Error resulting from the numerical 
determination of the bias level. 

Additional errors (e.g. count dependent bias) may arise 
because of a bad CCD or bad electronics, but dealing with 
those is beyond our scope. Fortunately we also had no 
significant dark current to deal with on any of the CCDs 
we have used so far. 

3.3. Generation of the error frame 

Our error frames are similar to a data-quality mask as 
used in many pipelines, but differ in that their values are 
not representative flags but actually numerical errors, with 
the exception of saturated and bad pixels which are rep- 
resented by simple flags (-1, 0). 

Each pixel {x, y) in an image / has an initial error 
5i{x,y) resulting from the photon noise, the bias noise 
(i.e. readout noise) plus the error in determination of the 
bias level: 



Siix,y) = 



< counts/ (a;, y) — bias/ 
gain/ 



bias/ 



(8) 
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where 

counts/ = flux of pixel {x, y) in image / in ADU, 
bias/ = bias of the image, 

photons r • r -L 

gam J — ^ ' = conversion factor, 

Cbias/ — we use the k{— 6)(T-clipped RMS of a 

suitable part of the overscan as an estimation 
for the bias noise (i.e. readout noise), 
'T'bias/ — number of pixels actually used for the bias 

determination. 
If using a small clipping factor (k < 3) for the determina- 
tion of the bias noise, (Tbiasj will be underestimated and 
therefore has to be corrected by a factor 1/Cc where 



K 

Cc = erf(K/V2) = ^ [. 

V 2vr J 



dk 



(9) 



Reproducible bias patterns can be determined with an ac- 
curacy only limited by the applied numerical precision, so 
their error may be neglected. 

We mark saturated pixels in the error frame by setting 
them to minus one, which will be dominant in any error 
propagation from now on to prevent the use of saturated 
pixels. 



3.7. Flatfield calibration 
3.7.1. Flatfielding basics 

In order to normalise the apparent photon sensitivity of 
all pixels in a single CCD frame, a calibration image (flat- 
field image) has to be built. In an ideal case this would be 
the image of an extended, homogeneous, flat, and white 
object at infinity. The apparent photon sensitivity results 
from geometric size, coating, and electronic properties of 
each single pixel, and, in addition, from the inherent prop- 
erties of the optics, and finally dust in the optics. Since 
there is no feasible way to get near the ideal case with 
dome flats (any additional optics which would project the 
nearby dome to infinity would again add features to the 
calibration image), we decided to stretch our efforts to 
improve sky fiats. The daylight sky is the ideal case for 
flatfield images, but is much to bright for broadband filter 
images. Therefore the suitable time for getting sky flats is 
restricted to dusk and dawn. The superiority of twilight 
flatfields over dome fiatfields is well known and e.g. dis- 
cussed in Bull (1991] ), and Mackay (1986 ). Nevertheless, 
dome flats and twilight fiats can be used effectively in 
combination as detailed in Sect. 3.7.3. 



3.4. Defective pixels 

We identify defective pixels by checking the specific CCD 
documentation as well as closely examining flatfields (for 
cold pixels, traps, and coating defects) and darks (for hot 
pixels) with a large spread of exposure times and counts 
if available. All sorts of defective pixels we set to zero and 
approximate them later (Sect. |3.9| ); the error is also set to 
zero. 



3.5. Photosensitive region 

Since the overscan parts of a frame are not needed any 
more, the frame is trimmed to its exposed part. We also 
truncate any bad, first or last rows or columns. On some 
CCDs the first row is broken and the edges of the photo- 
sensitive area behave differently from the rest: The effec- 
tive size of the edge pixels can be up to 20% larger than 
the size of an average pixel. 

3.6. Low counts pixels 

To ensure that only pixel values above an ignorance level 
will be used, we set low count pixels and their error to 
zero. This also accounts for the non linear range of a CCD 
due to a bad charge transfer efficiency (CTE) at low pixel 
charge^. Like the defective pixels (Sect. 3.4) those pixels 
will be estimated later (Sect. |3.9|), where it is possible. 



Since we do not have to deal with short exposures of bright 
objects in empty fields this procedure does not result in ig- 
noring the vast majority of pixels in a dataset. When deal- 
ing with empty, lo w count backgrounds we propose following 
e.g. 



3.7.2. Observation strategy for flatfields 

In order to minimise interfering effects like stars and pho- 
ton noise, wc need at least five fiatfield images per filter 
used, fulfilling these additional constraints: 

— The individual fiatfield images should have count lev- 
els as high or higher than the average of the observed 
object in a single frame, as their noise adds to the data 
(see Sect. 2.2 and 3.1). Nevertheless, fiatfields with less 
counts are preferred to having less than five flatflelds. 

— The exposures have to be long enough to avoid resid- 
uals of the shutter (shutter pattern|^, but still short 
enough to avoid saturation. If the shutter movements 
show a predictable time dependency, the fiatfields can 
be deconvolved from the two-dimensional shutter func- 
tion in a simple way as proposed by ^urma (1993 ). 
Unfortunately this is not the c ase with our obscrva - 
tions for the WeCAPP project (|Riffeser et al., 200l|) . 
Since the building of suitable flatfields already is the 
most (human) time consuming part of the standard 
reduction phase, we do not go into building shutter 
patterns for each individual flatfield, but simply ex- 
clude underexposed flats and try to manage with the 
remaining. 



McLean (1997) and adding a "fat" zero by preflashing the 



CCD. 



^ Assuming an average flatfield charge per pixel of 160 000 
electrons the shutter pattern will exceed one <j photon noise, if 
the exposure time is longer than 400 times the shutter move- 
ment time (opening plus closing) of a non photometric shutter. 
E.g. for an iris type shutter and a total shutter movement of 
10 ms the exposure time has to exceed 4 s just to have the ad- 
ditional shutter error not bigger than the photon noise. Since 
the shutter movement is a systematic effect, the combination 
of many short time flatfields will even enhance that error by in- 
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— Flatfields should be taken of a blank field to have less 
stars to be removed. 

— The telescope focus should be well adjusted to min- 
imise the number of pixels affected by stars. This is 
difficult for dusk flatfields, because one cannot make a 
focus series of a star before starting with the flat series 
for obvious reasons. 

— All flatfields of one series should have some small off- 
sets in orientation, so the Ku-clipped median of the 
series will be clean of stars. We found an offset of 30" 
to be sufficient]^. 

— Diffuse light contamination (i.e. reflections from inside 
the dome) has to be avoided. Neglecting this effect 
results in a mixture of dome and sky flat with improper 
illumination^. 



As one can easily estimate, using the Tyson fc Gal (1993 ) 
twilight formula, it is impossible to get five flatfield images 
(following above constraints) for more than one filter, if 
the time needed for CCD wipe and readout exceeds three 
minutes. For those cases we found the restriction to only 
take flat series for one filter per twilight, and to alternate 
the filters for each twilight period, to produce better re- 
sults than the combination of multiple insufficient flatfield 



3.7.3. Flatfield calculation and calibration 

The final fiatficld frame is built by combining individual 
flats following this standard recipe: 



1. We follow the steps described in Sect. 3.1 to 3.6 for the 
individual flatfield images to build appropriate fiats / 
and their correspondent error image Sj. 

2. Then we normalise all flats by dividing them with the 
median C/ of a central area of the CCD 



F = I/Ci 



(10) 



We determine the median only of the centre quarter of 
the full frame in order to minimise the effect of differen- 
tial illumination gradients when combining the flat^. 



creasing their fraction of all used flats and therefore amplifying 
their impact on the resulting flatfield. 

^ Our blank fields are clean of bright stars within the field of 
view of the cameras used (always less than 17 x 17 arcmin^). 
Since we also use only small, Im-class telescopes there is no 
considerable light contamination beyond the 30" limit of mod- 
erately bright stars or galaxies at average flatfield exposure 
times. 

^ One can check the impact of this light pollution effect by 
turning on a weak dome light while taking an image in a new 
moon night. We have improved considerably the situation at 
the Wendelstein telescope by painting the interior dome surface 
black. 

* The gradient due to vignetting within this central area is 
median level ± 3% in our images. A higher gradient may still be 
acceptable as long as it is guaranteed that the "normalisation 
median" lies within a well populated region of the normalisa- 
tion region's distribution. 



3. To get rid of stars, cosmics, and differential illumina- 
tion we build a KCT-clipped mean of the normalised flat 
frames F 



F{x,y) 



E F{x,y) 

F 



(11) 



^used 

For each pixel the outliers are excluded by Kcr-clipping, 
but we use the median instead of the mean as refer- 
ence for the selection procedure because the median 
is less sensible to occasional outliers; the mean of the 
remaining pixels yields the final calibration factor. We 
got the best results (least residuals) with a, hi — 1.0, 
which will preserve no more than two fiats for most 
pixels. With just five flatfields for combination a big- 
ger K left residuals at the 3% level. 

4. We control the result of step || by inspecting all con- 
trol frames i^controi — F/F. They should neither show 
any signal < 1, like holes or shadows around stars]^ 
nor illumination gradients. Both effects would indicate 
residuals still left in the median flatfield. 

5. Now we exclude those flats which cause residuals, and 
repeat from step ^ until there are no residuals left in 
the control images. If this would lead into having too 
few flats, we add the median flats of the days before 
and/or after to the median procedure. 

6. Since the pixel-to-pixel sensitivity variations are of- 
ten exclusively due to the CCD itself and therefore 
only a weak function of time, one can enhance a flat- 
field (i.e. improve its signal-to-noise ratio) by combin- 
ing the high spatial frequency, pixel-to-pixel variation 
of several twilight periods (or of very high signal-to- 
noise dome flats) with the smoothed flatfleld for the 
observed night, whose low spatial frequency informa- 
tion can change from night to night due to dust shad- 
ows etc. 

(a) We smooth each median flat pixel within a box 
smaller than the typical projected size of the short- 
time-scale variation due to dust: Each pixel (xq, j/o) 
is smoothed according to 

box 

Y.F{x,y) 

Fs{xo,yo) = — , where (12) 

^box 

nhoK = number of pixels in the smooth box. 

(b) The pure pixel-to-pixel variation is given by 
Fp{x,y)=F{x,y)/Fs{x,y) . (13) 

(c) The pixel-to-pixel variation flats are now combined 
by building a Kcr-clipped, median referenced mean 
as shown in step ^: 

"used _ 

E F^ix.y) 

Fp{x,y) = ^ . (14) 

^uscd 



^ Since the pollution of a flatfield image with stars or cosmics 
is always an additional signal, the origin of residuals in the 
control frames can be identified: Residuals with a signal > 1 
are correctly removed features of an individual flat. Residuals 
with a signal < 1, which also should be perceivable in multiple 
control frames, must originate from the median flatfield. 
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(d) The enhanced flatfield for the night can now be 

calculated via F{x, y) = Fs(x, y) Fp(x, y). 
This procedure will only be applied, if the observa- 
tional circumstances will lead to a gain in the signal- 
to-noise ratio of the resulting flatfield. A comparable 
technique has been used to build high quality HST 
flatfields (Ratnatunga et al., 1994). 
It is possible to check the gain (after a change of de- 
tector or an observational gap) using following corre- 
lation: 



1 



Ci 



'IIF, 



, where 



(15) 



F,Ci,F = see steps ^ and ||, 



'F/F 
^I/F 



RMS of all pixels of a control image F/F, 
RMS of an alternative control image I/F. 
This is only an approximation neglecting the readout 
noise and assuming identical pixel sensitivities. 

8. All data frames are divided by the finally accepted 
flatfield: i{x,y) = I{x,y) / F{x,y). 

We are still considering options to build a full automated 
flatfield evaluation procedure analysing the control im- 
ages. 

The whole flatfielding procedure can also be performed 
by standard astronomical data reduction software but 
without error propagation. 

3.7.4. Error propagation 

The statistical error is propagated as follows; all approxi- 
mations are only given to illustrate the impact of the re- 
spective reduction step and assume a normalised flatfield 
flux « 1 and negligible pixel-to-pixel and image-to-image 
variation of the error: 

1. Normalisation, error of pixel {x,y) in normalised flat 
F: 



6F{x,y) = F{x,y) 



' Sijx.y) 
. Iix,y) 



where 
F{x,y) 

I{x,y) 

Si{x,y) 

Ci 

5c. 



flux of pixel (x, y) in normalised flat F , 
flux of pixel (x, y) in bias corrected flat /, 
error of pixel (x, y) in bias corrected flat / 
defined in Eq. H, 



normalisation factor, see Sect. 3.7.3, step | 
error of normalisation factor, we use a 
standard error of the median 



\/»imcd 



nxncd = number of pixels used to build the median. 
An uniform CCD, homogeneously illuminated flat- 
fields, and a suitable (large) tailored area for the 



determination of C/ altogether will lead to a negligi- 
ble 5c I , which can be seen via 



SF{x,y) 
Sc. 



5i{x,y) 
Ci 

Si{x,y) 



1 



H , assummg 

^mcd 

and Ci ~ /(x, y) . 



Ci -y/J^mcd 

If using Kcr-clipping CTmod has to be corrected as shown 



in Sect. 3.3 



usmg Eq. ^. 

2. Building of median flatfield, error of pixel (x, y) in me- 
dian fiat F: 



F 



erf 



5f 



(17) 



^uscd ^total 
2 K>3 5f 



"total ^/"total 



, where 



^used 



remaining number of fiats after clipping 
used for a specific median clipped mean 
pixel, 

total number of fiats used for clipping, 
see Eq. |9[ and 
K, — clipping factor. 

Ignoring the effect of using preselected, Kcr-clipped pix- 
els to calculate the mean calibration pixel would lead 
to a significant misestimation of the resulting error for 
small n. The assumption of normally distributed values 



J^total 



is crude but still fair (see Sect. 6.1). 
3. The error of the fiatfield enhancement (if applied): 

(a) The error of a smoothed pixel (xq, j/o) built by av- 
eraging independent pixels is 



SfX^q^Vq) = 



/box 

li:5l{x,y) 



(18) 



"box 

ox 

(b) The error for the pixel-to-pixel fiatfield is given by 



/ 6p{x,y)y 1 


( ^F, 


y)\ 


\F{x,y)) 


\Fs{x, 


y)J 



(19) 



(c) Combining the pixel-to-pixel fiats yields (as in 
step H) 




E SI (x,y) 



erf (^^) V^Wd 



(20) 



"total 
K>3 Sf„ 



"total V"total 

(d) The error of the enhanced fiatfield can be calcu- 
lated with 

Sp{x,y) = F(x,2/) (21) 
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K>3 

dp \l 

n-hox '^total 

4. Flatfield division, error of pixel {x,y) in flatfield cali- 
brated image /: 



The first approach does not work at all if there are no 
multiple images available or the sample of images to be 
stacked has different observational features (variable sky, 
extinction or PSF). Aligning images will always spread 
and diffuse cosmics and therefore obscure them. The latter 
technique gets into trouble with noisy images, undersam- 
pled images and multiple- pixel cosmics for obvious rea- 
sons. 



Trainable cosmic classifiers i.e. as described in salzberg 
et al. (1995| ) have the advantage of also being applica- 



F{x,y) 



(22) 



Minor Bystematic errors are neglected here: 



— The flatfield response of a CCD is a strong function 
of colour. This results in a systematic error when cal- 
ibrating stars with colours different from the sky on 
an image, and gets worst when the compared objects 
have highly different colours. 

— The geometric distortion introduced with the variation 
of CCD pixel size is ignored in our flatfielding proce- 
dure. A position estimate will have a systematic error 
according to this, if determining positions of objects in 
undersampled images (e.g. due to extraordinary good 
observing conditions and therefore very sharp PSFs). 

— Since we ignore the individual geometric sizes of CCD 
pixels the integrated photometry of a flatfielded image 
may be corrupted (not exceeding 0.5% in a single pixel 
of our images). 

All those errors could be compensated, but will have at 
most a minor (but detected) influence on our data because 
of OIS, relative profile fitting photometry, dithered image 
stacks and error propagation (Sect. 5^, 5^, and 5.5). 



3.8. Filtering of cosmic rays 

There are two major reasons, why we have to correct for 
the image contamination by particle events (so called cos- 
mics): We use small, Im-class telescopes for our obser- 
vational projects (e.g. Riffeser et al., 2001) and therefore 



have integration times of half an hour or longer. We want 
to find variable sources and establish light curves for every 
pixel of an observed field. So we have to identify individual 
cosmics in single frames automatically at least, or better 
clean the images of cosmics and account for the error this 
procedure might introduce in a stack of images. 

3.8.1. Common and literature filters 

Existing filtering techniques for particle events either rely 
on median stacking of multiple well aligned images (i.e. 
see Windhorst et al., 1994) or compare each pixel value 



blc to undersampled data but rely on subjectively de 
fined training sets which are difficult to create for a large 
spread of different telescope, camera, and detector config- 
urations in addition to the great range of observing condi- 
tions. An interesting new approach is presented in Rhoadg 
(2000D . Since this method relies on an accurate PSF and 



sky determination and the author does not refer to possi- 
ble problems in heavily crowded fields we did not use it. 
Furthermore, this technique is in principle only sensitive 
to single pixel events, which we found to be not the com- 
mon case. In fact most cosmics seem to have a major-to- 
minor axis ratio|^ greater than two. Multiple-pixel events, 
which can be filtered with our technique (see below) in one 
pass, have to be detected iteratively and with decreasing 
efficiency with the Rhoads technique. 



3.8.2. Gaussian filter 

We apply a straightforward Gaussian filter to every single 
image: We fit five-parameter Gaussians to all local max- 
ima of an image. If the width along one axis of the fitting 
function is smaller than a threshold (which has to be cho- 
sen according to the PSF) and, in addition, the amplitude 
of the fitting function exceeds the expected noise by a fac- 
tor (which has to be chosen according to the additional 
noise due to crowding, see Sect. 3.2 for details), we re- 
place the pixels with the fitted surface constant, where 
the fitting function exceeds this constant by more than 
two times the expected photon noise. In the following we 
describe the algorithm in detail: 

1. Because of code speed improvements which rely on 
some symmetries in the fitting function the tested cos- 
mic candidate has to be in the centre of a 7 x 7 pixels 
array. In order to be applicable also on the first and last 
three rows and columns we add a border surrounding 
the exposed frame filled with zero value pixels. Since 
we do not want to lose too much of the images when 
shifting them later (Sect. ^ we chose to enlarge the im- 
ages not only with a three pixels but with a 20 pixels 
border instead. 

2. Now, first we search for all local maxima {xo,yo) in 
the image, but ignore those with either a large erroip] 



with the median of its neighbours and define pixels with 
a sharpness ratio above a deliberately set value as cosmic. 



We have checked this with the control output of our filter 
code. It gives the major and minor axis full width half maxi- 
mum of the Gaussian fit function for every cosmic replaced. 

We found 7 s:; 3 to be an empirically suitable factor. 
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(^6'^{xo,yo) > 7^ ^'^"tahi''^"^ ) more than two satu- 
rated neighbours'^ or with less than four (of eight pos- 
sible) valid neighbours]^. 

Then we perform a propagated error weighted, 
least-squares fit, assuming a five-parameter, two- 
dimensional Gaussian fitting function centred on these 
local maxima (xq, yo) coordinates in 7 x 7 pixels subar- 
rays: We determine a surface constant C, an amplitude 
A, a rotation angle a, a major and a minor axis full 
width half maximum (xt^hm a-nd ?/fwhm) of the fitting 
function /gauss giving the flux of a pixel (x, y) 



/gauss (a^j^) — C + 



(23) 



Aexp 



-4 In 2 



fwhm 



^fwhrn 



where 



x'{x,y) — (a; — xo) cosa -f (y — yo) sina, 
y'{x,y) = (y - yo) cosa - (x - zo) sina. 

All Gaussians with an amplitude of ^umit times the 
propagated error of the centre pixel and a full width 
half maximum in one axis smaller than a limiting 
FWHMcosmic are defined as cosmic: 
{A > iiimit ^(a;o,yo)) A 

[(a^fwhm < FWHMcosmic) V (yfwhm < FWHMcosmic)] 

We have to perform a sanity check on the fitting func- 
tion: The surface constant C and the amplitude A must 
be positive; of the fit must be close to unityf^: 



C > A A>0 A Xfit - 1) where 



(24) 



Xfit 



the reduced of the fit 



fitbox 

E 



'^dof 



«dof 

degrees of freedom for the fit, 
(7 X 7) - 5 here. 



(25) 
(26) 



6. Now we substitute pixels (x, y) where the Gaussian 
fitting function /gauss gives a value larger than the as- 
sumed signal plus two times the assumed photon noise 

f gauss (x, y) > C + 2^ 

7. The new error of the substituted pixels is set to 



l6j{x,y)+xlt 



C 
gain 



where 



(27) 



Cosmics candidates close to saturated pixels need a special 
treatment, see step ^. Because we always consider the possibil- 
ity of blooming (see Sect. 3.1), only for pixels with at least three 
saturated marked pixels (of the eight possible neighbours) one 
pixel (of the four directly adjacent) is really saturated. 

The Gaussian fit of step ^ gets unstable when fewer than 
half of the pixels adjacent can be used in the fitting algorithm. 

This corresponds to a compatible fitting function and cor- 
rect error weights. 



5j{x, y) = propagated old error of pixel, 

Xfit = accuracy of the Gaussian fitfunction 

(x ~ 1 for a perfect fit). 
This will be large enough to prevent an incautious use 
of the replaced pixels in the following. 
We then repeat this procedure for the areas with cos- 
mics found beginning with step ^ until no more cosmics 
are found. 

Finally we try to find and replace cosmics near satu- 
rated pixels with a similar, just in some details more 
sophisticated technique: The Gaussian fit starts cen- 
tred on the saturated region but the centre position 
is added to the list of free parameters. We ignore the 
saturated region for the fitting and the replacement 
procedure. For overall stability reasons we have to use 
closer constraints for the sanity check. Since saturated 
pixels due to cosmics will not be treated at all, because 
they are flagged as "dominant bad pixels" , step ^ might 
be readjusted, if dealing with shallow-well CCDs; sat- 
urated regions can be replaced, but then saturated ob- 
jects may be mistaken for a cosmic. 



We found that a iumit = 8.0 and a FWHMcosmic = 1-5 
works fine in any well sampled image. However, in some 
extraordinary good seeing images (with an average stel- 
lar PSF FWHM < 2.0 pixels) we had to specify a Hmit- 
ing FWHMcosmic = 1.3 to avoid the deletion of stars. All 
fixed fitting and substitution constants were adjusted in 
order to get an accurate and reliable filter for cosmic rays 
for all our images. The sensitivity parameters inmit and 
FWHMcosmic have nevertheless to be adjusted in general 
to the observational and object properties to reach the 
best compromise between false alarm and fail detection 
rate (Sect. O). 



3.9. Approximation of bad pixel areas 

Pixels with value and error set to zero will now be re- 
placed. We use a distance and error weighted linear ap- 
proximation of the closest neighbours. The fitting box is 
selected as small as possible with the restriction that more 
than 2/3 of the fitting box pixels minus the centre pixel 
must be valid pixels and the fit box may not be larger 
than an arbitrary limit which we set according to the spa- 
tial resolution of a specific imaging system. If even the 
largest possible box does not apply to the first criterion 
the pixel is considered as isolated and not replaced at this 
point^. Each replaced pixel (xo,yo) of an image / gets 



Isolated pixels still can be replaced with the mean value 
of corresponding pixels in the remaining images of a stack 
(Sect. ^). We do not apply this method, which relies on a per- 
fect photometric alignment of the to-be-stacked frames before 
any image convolution has been applied. If we need a measure- 
ment in the lost area, we build a separate stack without the 
image with isolated pixels in the interesting region. 
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an error calculated from the individual pixel errors of the 
fitting box 



"used 

^uscd 



where 



(28) 



'^used = number of pixels used for fit, and 
T^fj^ = defined according to Eq. |2^. 

A larger fitting box has fewer close base points and 
therefore raises the uncertainty of the linear approxi- 
mated substitute, so we use an average error of the input 
parameters times the quality of the fit (xfit) and not 
only th e er ror of the calculated value (cx Like 
in Sect. 3^ , step this prevents an incautious use of the 
replaced pixels in the following. 



4. Position alignment and stacking 

Up to now our reduction pipeline can be applied to any 
image regardless of its scientific application. When shift- 
ing images in order to stack them the PSF and even the 
flux may be altered. It is important to conserve the PSF, 
because the image convolution approach that we adopt 
performing differential photometry relies crucially on it. 

In any case the alignment of images follows a four step 
procedure: 

1. First we determine the coordinates of reference objects 
in every image, 

2. then we calculate the coordinate transformation to 
project an image onto the reference frame, 

3. subsequently we project the images into the reference 
frame coordinate grid, and 

4. finally we stack the images. 

4.1. Position of reference stars by PSF-fitting 

In order to obtain the coordinates of reference objects 
in all images we perform interactively a seven-parameter 
Gaussian fit. We begin with the reference frame: About 20 
stars with a high signal-to-noise ratio and well distributed 
over the frame would be sufficient but we use 50, because 
stars close to the frame border may be missing on some 
images. We continue with selecting at least one reference 
object in every image manually, the rest will be found au- 
tomatically]^. The lists with the reference objects in the 
reference frame and the first reference object of the un- 
shifted images are used to recognise the reference objects 
in each image, to determine their position, and finally to 
calculate the projection parameters. 



In case the imaging device provides an accurate World 
Coordinate System (WCS) information all of the reference 
stars could be found automatically. 



4.2. Translation of image coordinates 
of a linear projection 



determination 



With the telescopes and cameras used in our observing 
campaigns, we found a linear relation to be sufficient. We 
easily match 50 stars all over a 8' x 8' field within 1/20 " 
in the mean. Since there was no significant optical field 
distortion, it was not necessary to use a non linear relation. 
We determine a 2 x 2 linear matrix and a two-dimensional 
translation vector 



y 



ail ai2 

121 122 



(29) 



with a least-squares fit. It matches the positions of refer- 
ence stars in the reference system with the positions in 
the unshifted image with this six-parameter relation. 

4.3. Shifting of images 

The tech nique of Variable-Pixel Linear Reconstruction or 
drizzling (Fruchter fc Hook, 1997| ; [Hook fc Fruchter, 1997 ; 



Mutchler fc Fruchter, 1997 ) offers the possibility to add 



images while preserving both the flux and the PSF. In un- 
dersampled images one might even enhance the resolution 
and therefore gain signal to noise. Unfortunately this tech- 
nique requires some image properties to be applicable and 
these are very difficult to obtain with ground-based tele- 
scopes: There must be no variations in sky, extinction and 
PSF, and there should be a uniform spatial sampling in the 
sub-pixel pointing. If those requirements are missed, the 
flux will still be preserved, but the PSF may get very dis- 
torted. So aperture photometry might still work very well, 
but since we use PSF convolution and a PSF-dependent 
photometry, we had to think for an alternative with less 
observational constraints. 

Our shift algorithm preserves PSF in unstacked and 
sometimes undersampled frames. We found that a 16- 
parameter, 3'''^-order polynomial interpolation with 16 
pixel base points does apply to our needs. A 2"'*-order 
polynomial still smoothes the images, whereas a 4*'^-order 
polynomial does no better PSF conservation compared to 
the 3"^-order polynomial. Since the number of parameters 
of the polynomial is matched with the input base points, 
no least-squares fit is needed; the polynomial can be cal- 
culated analytically. 

The flux interpolation for non-integer-value coordi- 
nates {x, y) is calculated with a polynomial 



3 3 



P{^^ y)=^^ "■ij^^'y' ' where 



(30) 



i,i ~ index for subscript, and exponent for superscript. 
For a region with 4x4 pixels this yields 16 linear 
equations 



p{xk,yk) = Ii^k,yk) , where 1 < fc < 16, 



(31) 



so the coefficients — aij{I{xk,yk)) can be calculated 
by solving the matrix equation. The error is calculated 
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using Gaussian error propagation 
6p{x,y) = 



3 3 



\ 1=0 j=0 



(32) 



dai.j[i{xk,Vk)) 



4.4. Stacking images 

In order to stack images we simply sum up the aligned 
frames, but exclude outliers (satellite polluted, occasional 
bad seeing or bad focus images) and if convenient split 
into groups (i.e. if many observations available of one night 
or very variable PSFs etc.). When searching for variable 
sources one always has the choice of trading time reso- 
lution for a deeper magnitude limit. One could also get 
both, if looking for periodic variables, but then one has 
to test every possible period with different stacks and one 
may also have to stack images with very different obser- 
vational properties. Pixel defects still marked with zero in 
any frame (like saturated or not interpolatible pixels) are 
set to zero in a stack. 

The error of a pixel (x, y) in an image stack S consist- 
ing of rij images propagates as 



Ss{x,y) = 



(33) 



5. Detection of variable objects 



In order to detect variable sources we are using DIA 
[Difference Image Analysis)^ a method proposed by 
^omaney fc Crotts (19961 ). The idea of DIA is to sub- 
tract two positionally and photometrically aligned frames 
which are identical except for variable sources. The result- 
ing difference image should then be a flat noise frame, in 
which only the variable point sources are visible. 

A crucial point of this technique, apart from posi- 
tion registration (Sect. ^) is the requirement of a perfect 
matching of the PSFs between the two frames. For this 
purpose we apply a method, called OIS [Optimal Image 
Subtraction) developed by Alard fc Lupton (199S| ). 



5.1. Photometric Alignment 

Although OIS implements the photometric alignment 
from the reference frame to each other frame, it is useful 
to align all frames photometrically before matching the 
PSF. This ensures that the whole light curves are photo- 
metrically calibrated to a standard flux. 

All frames S[x,y) are connected to a flux standard 
frame So[x,y) by[^ 

So[x,y)^a S[x,y)+b = I[x,y) . (34) 



To be consistent with the Alard & Lupton (1998) notation 



we now call images again plain /. 



The scaling factor for different exposure times and atmo- 
spheric extinction a and background sky light b are deter- 
mined in a simple and crude way. As a first step we remove 
all obvious bright stars from our field and replace them by 
a plane representing the surrounding background level. By 
replacing in both frames each pixel with the median count 
rates of 21 x 21 pixels subsections, the influence of a dif- 
ferent PSF between the frames on the determination of a 
and b is minimised. With these replaced pixel images we 
calculate values for a and b by solving the least-squares 
problem. 

The error frame is calculated using Gaussian error 
propagation 



Si[x,y) = ^ S'2(x, y) 51 + 5|(x, y) + 5\ 



(35) 



5.2. Convolution with differential background 
subtraction 

One advantage of OIS is the possibility to fit differential 
background variations simultaneously with the PSF be- 
tween the frames. 

Including a background term the convolution equation, 
which transforms the PSF with smaller FWHM of the 
reference frame R to the PSF of an image /, is of the form 

I[x, y) « R[u, v) (g) K[u, v) -f bg[x, y) = R[x, y) , (36) 

where [R ® K)[x, y) — J2 + 2/ + v)K[u, v) . 

u,v 

The convolution kernel K[u, v) and the background term 
bg[x,y) are decomposed into basis functions 

n 

K[u,v) — aiBi[u, v) , and 

n+Tihg 

bg[x,y) = E a.xP'y'^' , 

i—n+l 

where n is the total number of coefficients of K[u, v) and 
ribg is the corresponding number for the background term 
bg[x, y). The exponents pi and qi are fixed integers. 

K[u, v) is determined by solving the least-squares 
problem: 

= E^[(i?§5if)(a;,y) + 65(^,y)-^(2;,2;)]' (37) 



By setting = these equations transform into 

E"'E^C'.(a:,2;)C,(a:,y) (38) 

i x,y ^-y 

= ^^-I[x,y)Cj[x,y) , where 

x.y ^^y 

R(u^ v) (8) Bi(u^ v) z = 1, . . . , n 
C,[x,y) = { xP^yi' i^n + 1,..., (39) 

n + Ubg 
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The problem is reduced to the solution of the following 
matrix equation for the Oi coefficients 



After the convolution the difference frame D is com- 
puted by subtracting the R frame from the I frame 



at My ^Vj , Ma = V (40) D = I-R 

i 

where the matrix elements are defined according to 

= E ^C^i^'V) ' (41) 

X.y 

= E:;^^(^'y)^^(^'y)- (42) 



(45) 

The calculation of the error frame is done according to 
Gaussian error propagation 



Siiix, y) = /E ^Iv ^l+u,y+v and 



x,y 



^oix, y) = J SUx, y) + J|(x, y) . 



(46) 
(47) 



5.3. Application to data 

We adopt Gaussians modified with polynomials of order p 



as a suitable kernel model as proposed by Alard & Lupton 

K{u^ v) — a.iBi{u^ V 

i 

=E 



^^+1,^ VI Pi -3 



E E "'j'^ "^^'^ 



; j=o k=o 

In order to get an optimal kernel, which can be of a compli- 
cated shape, we use the combination of three Gaussians 
( Alard fc Lupton, 1998| ), each with a different width a. 
This leads to the following 49 parameter decomposition of 
the polynomials of K{u,v): 



(Tl = 1, pi = 6 
= 3,P2 = 4 
0-3 = 9,_P3 = 2 



(ai 4 

(a29 ■ 
(044 ■ 



028"^) 

+ a43U^ 

+ a49M^ 



(43) 



Additionally ribg = 3 parameters are implemented to fit 
the background 



bg{x, y) = flso + a^ix + 052?/ 



(44) 



To get rid of the problem of a spatially var ying PSF 
over the whole a rea of the CCD, as described in Tomaney 
fc Crotts (1996|) , we divide the images in regions, each 



containing 141 x 141 pixels. In each of these regions a 
locally valid convolution kernel is calculated. Since K(u, v) 
has a box size of 21 x 21 pixels and the regions to be 
convolved are overlapping, each region comprises 161 x 
161 pixels. For all bad pixels (marked as 0) the convolution 
is not done, these pixels remain 0. 

The calculation of the matrix Mij is the most time 
consuming part of the convolution. The matrix Mij of the 
reference frame R is calculated once and can be used for 
all images. To enable this timesaving approach we take 
the error Ux.y which enters the calculations always from 
the error frame of R {(Jx.y = ^Rix^y)). Therefore only 
the calculation of the vector Vj has to be done for each 
image/reference frame pair. 

Bad pixels in the frame / would lead to an error be- 
cause they are not marked in the frame R. To compensate 
this the calculation of the matrix is redone for these pixels 
and then subtracted from the original matrix. 



5.4. Source detection 

We developed a standard star finding algorithm to detect 
sources in the difference images: We smooth the image by 
replacing each pixel with the mean of five pixels of a cross 
shaped area and tag all local maxima. Then we fit a simpli- 
fied Moffat function (|Moffat, 1969| ) A[l + s(a;2 + y^)] '^'^ 
to these local maxima in the unsmoothed image. After ex- 
cluding all bad fits we fit a rotated Moffat function to the 
remaining maxima: 



/moffat(a;, y) ^ A[l + {sxx') + {syy') 



(48) 



where 

x' = {x — xq) cos a + (y — yo) sin a, 
y' = (y - 2/o)cosa - (a;- a;o)sina. 

a denotes the rotation angle, A the amplitude and the 
pair xo,yo the central coordinates of a stellar PSF. The 
rise of the wings of the PSF is given by the parameter /3, 
whereas s, Sx, and Sy specify the width of the function. 
We include the errors taken from the error frame in the 
nonlinear least squares fit considering the error frame to 
weight the count rates obtained in the frames. 

Minimum and maximum expected FWHM of the PSF 
and minimum and maximum of P have to be chosen ac- 
cording to observational conditions. To distinguish be- 
tween noise and real sources a threshold factor t, is in- 
troduced; t gives the ratio of the parameter A and the 
background noise. All sources below a certain threshold 
(i.e. t = 5 for the WeCAPP project) are regarded as noise. 
Because difference images can comprise negative sources 
the images are inverted after one detection cycle. The 
whole detection procedure is then redone on this inverted 
frame. 



5.5. Photometry 

Photometry of the detected sources is performed with a 
profile fitting technique. We choose reference stars in the 
CCD field to obtain the information about the PSF of 
any particular frame. These stars should be bright and 
isolated enough to allow an accurate determination of the 
PSF. On the other hand they have to be unsaturated in 
any of the images. Since all frames are already photomet- 
rically aligned (Sect. 



5.1), the amplitude of standard stars 



is not used for flux calibration, so standard stars may even 
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be variable. We apply a Moffat fit (Sect. ^ Eq. ^) to 
some standard stars selected in such a way in each image. 
Keeping the slope parameters /3, s^, Sy, and the angle a 
of the best fit star only the amplitude A has to be deter- 
mined for each detected source in the difference frame by 
a linear least-squares fit using propagated errors. We inte- 
grate the count rates over the area of the now fully known 
analytical function of the PSF of the source to determine 
its flux. This minimises the contamination from neighbour 
sources. 

The fitting error of the amplitude 



1 



/mottat(a;,y) 
A S'i,(x,y) 



(49) 



is transformed into the error of the flux determination by 
multiplying 6a with the flux corresponding to A = 1. To 
account for the accuracy of the PSF-parameter determina- 
tion we multiply the resulting error with the Xstandard star 
(Eq. ^) of the standard star fit: 

^Flux — 5a Flux^^i Xstandard star • (50) 

6. Simulation tests with errors and error 
propagation 

To show the influences of a propagated error estima- 
tion we have performed image reduction with and with- 
out error propagation on simulated images, which include 



all the image properties of Sect. 2.1. These simulations 
were done in addition to the empiric examination of real 
data to have well defined input parameters and there- 
fore be able to decide whether the expected performance 
can be achieved. All test images represent closely one of 
our detector configurations: Ik x Ik CCD, 16 Bit ADC, 
gain = 3, bias = 200 ± 3 ADU, saturation = 65 000 ADU. 

6.1. Errors in flatfields 

We show the impact of error propagation by flatfielding 
artificial skylight images with both flatfields and images 
representing different flux levels. Tab. || gives the mea- 
sured relative errors i.e. the normalised standard devia- 
tion of the image. The realistic ~ 30k counts flatflelds 
case compared to the naive Si case (Eq. |^ without any 
further correction for flatfielding) gives an underestima- 
tion of 5% to 30% compared to the true errors. For the 
very low counts (but clean of stars and cosmics) case the 
underestimation is even 50%. The mean propagated error 



estimate, calculated as shown in Sect. 3.7, always differs 
less than 2% from the measured error. 



6.2. Errors with Gaussian filter for cosmic rays 

Empirical tests on real images were done by comparing 
the number of cosmics found in exposed images with that 
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Table 3. Normalised standard deviation of flatfielded 
artificial skylight images [%]: 6i denotes a perfect flat, 
with noise exclusively induced by the skylight image, i.e. 
the naive error; rows 10 to 60 in respect to the 5i row 
illustrate the impact on the error budget for different flux 
levels of the median of five fiatfield calibration image; ~ 
30 shows the realistic case of five artificial dithered flats 
comprising stars and cosmics. 



found in dark frames and visually examining both the un- 
filtered image and the difference of the unfiltered and the 
filtered image (e.g. for effects on bright stars etc.). 

In addition we have tested the reliability of our 
Gaussian filter with five simulated cases (Tab. |): A pure 
skylight image, a simple field with 500 plus (A'sat =) 
10 saturated and bloomed stars and two different sky 
levels, a crowded field with 100 000 stars, and a highly 
crowded field with 200 million stars; the positions of stars 
and cosmics as well as the flux and the orientation of 
cosmics follow uniform deviates; the flux of stars fol- 
lows a exponential deviate, which is a sufficient match 
to the luminosity function of our fields. Stars have a PSF 
FWHM « 2.6 pixel. The images are processed with our 
standard reduction pipeline (Sect. H) using a median- 
of-five simulated fiatfield with an average flux level of 
30 000 ADU per flatfield and the filter parameters of 
Sect. 3.S unless stated otherwise (detection thresholds 

tlimit = 8, FWHMcosmic = 1-5). 



We determine the false alarm rate by filtering the clean 
test images without any cosmics (Tab. I): For lO*' pixels 
with about 30 000 to 100 000 local maxima 0.3 to 1 x 10^ 
tests are performed. So the false alarm rate is given by the 
ratio of false occurances to number of tests. To determine 
the detection rate we put 500 artificial cosmics with fiat 
deviates in space, form, and energy into the test images 
(Tab. 1^): We identify and count the cosmics found. To get 
an accurate estimate of the performance these numbers 
still have to be compared with the photon noise, the noise 
induced by the object density, and the filter parameters. 

The highest false alarm rate is achieved with the 
crowded field. Here the total pixel-to-pixel variation of 
the image (photon noise plus objects' signal) exceeds the 
pure photon noise by a factor of 25. In the highly crowded 
field this excess is only a factor of 1 1 and can be compen- 
sated by setting tumit — 10. The false alarm cx A"sat in the 
simple, high sky field is due to the unawareness of satu- 
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Table 4. Parameters of the test images for the Gaussian 
filter for five different test configurations (fields) showing 
the level of background sky [ADU], numbers and maxi- 
mum fluxes [ADU] of stars and cosmics. 
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Table 5. Performance of Gaussian filter for cosmic ray 
events for the test configurations of Tab. |[ The false alarm 
and the detection rates with and without error propaga- 
tion frames; f indicates that the detection rate matches 
the expected rate according to iumit within 2 cosmics — 
0.4%. (a: tu^it = 10.) 



rational because of missing saturation tags without error 
propagation. Our tests show a false alarm rate < 0.2% in 
any (error propagated) case. We found that our fields are 
resembled closer by the smooth simple-and-high-sky field 
than by the crowded and even the highly crowded test 
fields. Therefore we are sure the false alarm rate does not 
exceed 0.01% in our real images. However, a false alarm 
resulting in a deletion of a true source will not lead to 
wrong photometry because of the high error assigned to 
replaced pixels (Sect. |3.8| ). Under certain circumstances 
(bad sampling, small stacks) it may lead to large error 
bars, but the result is still reliable within those. 

The expected detection rate (97% to 99%) is achieved 
nearly with all test images including error propagation. 
Even without error propagation the detection rate is still 
very close to the expected one. The worst case is again the 
(not highly) crowded field where the object induced noise 
exceeds the photon noise by far. Here the expected rate is 
missed by 0.008 = 4 (of 500) cosmics. 
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Fig. 1. Accuracy of the projection: position differences 
(Ax — x[ — X2, Ay = 1/1 — 2/2) of 70 stars after projecting 
the coordinates of one frame {xi,yi) to the other frame 
(2^2,2/2). 



6.3. Accuracy of the linear projection 

The projection (Sect 



4.2) was tested with two simu- 
lated, not perfectly aligned (shifted, rotated and rescaled) 
frames. It was calculated to match the position of 70 bright 
stars in these frames. The position differences are always 
below 0.05 pixels (Fig. |l]). This refiects the principle ac- 
curacy limit^Jdue to the size of the corresponding fit box 
(20 X 20 pixels). 



6.4. Errors visible after alignment - a snapshot 

To give an impression of error features which would be 
neglected by just considering the cleaned image, but still 
visible in a propagated error image, we present image and 
error image of one hour light on a small telescope of the 
dwarf galaxy EGB 0427 -f 63 (Fig. |). Despite the fact that 
the images were dithered there are still features of the 
flatfield (dust rings) clearly visible as well as the impact 
of CCD defects and a huge amount of cosmics. 



6.5. Errors in convolution 

We tested the accuracy of the convolution with 19 pairs 
of simulated frames: The reference frames with a FWHM 



As stated in Sect, p.q saturated pixels need a special treat- 



ment 

19 



Spatially extremely undersampled images, leading to peak- 
shaped PSFs, still can restrict this principle accuracy limit to 
one pixel, but this was not a required test case. 
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Fig. 2. The dwarf galaxy EGB 0427 +63: One hour Ught on 0.8 m WST telescope, best seeing and lowest sky 20 images 
of four days stacked; left: image; right: error image. 



of the stellar PSF of 2.4 and five times more flux than 
the comparison frames with a FWHM of 3.0, all with over 
100 000 stars and different background levels. According 
to Sect. 5.2 the reference frames are convolved to match 



the PSF and the background level of the comparison 
frames. The convolved frames are subtracted from the 
comparison frames and the result divided by the expected 
RMS errors, as derived from error propagation. This gives 
the ratio of expected photon noise and measured noise. 
The histogram of such a ratio frame matches a Gaussian 
with (7 = 1 almost perfectly, which means that the ex- 
pected photon noise fits the measured noise. This shows 
that the OIS method can be applied to very crowded fields 
like M 31 and gives residual errors at the photon noise level 
(Fig. I). 

6.6. Errors in interpolations - PSF-fitting photometry 

We performed PSF-fitting photometry in simulated im- 
ages as described in Sect. 5.5 using a Moffat fitting func- 
tion (Eq. HI) on bright stars, and (together with OIS, us- 
ing a high signal-to-noise reference image) on faint variable 
sources in a crowded and a highly crowded field. We found 
that per pixel propagated errors compared to estimated er- 
rors greatly enhanced the reliability of any fit. Extensive 
testing with simulated images comprising different ob- 
servational features shows that this is especially due to 
the treatment of defective pixels, cosmics and saturation 
(pixel defects) . If we want to avoid the labour of full error 
propagation we nevertheless have to use masks to get rid of 
these pixel defects. But then long time series will diminish 
our field, because the defects will spread (due to dithering, 
minor misspointing, different detectors, random position 
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Fig. 3. Histogram of the pixel values of a simulated differ- 
ence image divided by the expected RMS errors. The solid 
curve is a Gaussian with a — 1. We calculate the reduced 
chi-square of 19 different simulated images in the range 
between -3 and 3. The median is 1.1, which means that 
expected and measured errors match almost perfectly and 
that the residuals in the OIS are at the photon noise level. 



defects etc.). Furthermore, the (Eq. of a fit has a 
valid meaning only for fully propagated errors: « 1 im- 
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plies a correct measurement within purely noise induced 
errors and correctly flagged pixel defects; > I indicates 
systematic errors beyond noise and corrected pixel defects 
like blending of variable sources, missed or wrongly treated 
pixel defects. Renouncing full error propagation will shift 
e.g. bad flatfields from the recognised, high noise category 
to the undiscovered systematics regime. The calculated er- 
ror will not change, because we consider for the final 
error budget, but we would lose the chance to investigate 
those cases. If one cannot prove the origin of the uncer- 
tainties in a measurement, one can neither be sure of the 
measurement itself nor of the postulated accuracy. 

7. Summary 

We have presented a standard image reduction pipeline 
applicable to any CCD images, but with the intention to 
execute DIA, i.e. image convolution and relative photo- 
metry. We put a great emphasis on the importance of per 
pixel propagated errors. The necessity of building good 
flatfield calibration images and a way to obtain them 
was shown. We also presented a robust filtering technique 
for cosmic rays applicable to single, not too undersam- 
pled images. We discussed all image reduction issues of 
finding variable objects and measuring their variations 
in highly crowded fields: PSF matching by convolution 
of a reference image using the Alard fc Lupton (1998 ) 
technique, construction of difference images, detection al- 
gorithm for variable sources, and relative photometry of 
variable sources by a profile fitting technique. 

A paper about the application of this image reduction 
pipeline on a massive imaging cam paign of a part of th e 
M31 Bulge is accepted by A&A ( [Riffeser et al., 20011 ). 
Parts of the reduction pipeline have also been successfully 
applied to MUNICS data ( cosmic-filtering of MOSCA 
spectra and CAFOS images; [Prory ct al., 200 1[ ) and VLT 
FORS data (cosmic-filtering and image alignment of re- 
vised FDF frames, A. Gabasch, priv. comm.; image align- 
ment and OIS in the centre part of NGC 4697 using a dif- 
ference image built of narrow on-band and off-band line 
images, Mendez et al., 2001). 
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